Root_Atlas.rds (or get it by running through notebook 2, 3, 4, 5, 6, 7, 8, 9 & 10)
features.tsv.gz
Curated_Markers.csv
rm(list=ls())
# Set the working directory to where folders named after the samples are located.
# The folder contains spliced.mtx, unspliced.mtx, barcodes and gene id files, and json files produced by scKB that documents the sequencing stats.
setwd("/scratch/AG_Ohler/CheWei/scKB")
library(tidyverse)
library(Seurat)
library(RColorBrewer)
library(future)
#for 200gb ram
options(future.globals.maxSize = 200000 * 1024^2)
# Load atlas
rc.integrated <- readRDS("./Root_Atlas.rds")
# Simple QC label
rc.integrated$celltype.anno <- gsub("Putative Quiescent Center", "Quiescent Center", rc.integrated$celltype.anno, ignore.case = FALSE, perl = FALSE,
fixed = T, useBytes = FALSE)
order <- c("Quiescent Center", "Stem Cell Niche", "Columella", "Lateral Root Cap", "Atrichoblast", "Trichoblast", "Cortex", "Endodermis", "Pericycle", "Phloem", "Xylem", "Procambium", "Unknown")
palette <- c("#9400d3", "#DCD0FF", "#5AB953", "#BFEF45", "#008080", "#21B6A8", "#82B6FF", "#0000FF","#FF9900","#E6194B", "#9A6324", "#FFE119","#EEEEEE")
rc.integrated$celltype.anno <- factor(rc.integrated$celltype.anno, levels = order[sort(match(unique(rc.integrated$celltype.anno),order))])
color <- palette[sort(match(unique(rc.integrated$celltype.anno),order))]
# Plot atlas with celltype annotation
options(repr.plot.width=8, repr.plot.height=6)
DimPlot(rc.integrated, reduction = "umap", group.by = "celltype.anno", cols=color)
# Plot atlas with time zone annotation
options(repr.plot.width=8, repr.plot.height=6)
time_plt <- DimPlot(rc.integrated,
group.by = "time.anno",
order = c("Maturation","Elongation","Meristem"),
cols = c("#DCEDC8", "#42B3D5", "#1A237E"))
time_plt
# Set the identity of interest to cell type
Idents(rc.integrated) <- "celltype.anno"
# Set the expression matrix used for DE to batch-corrected and scaled one
DefaultAssay(rc.integrated) <- "integrated"
# run differentially expression (DE) analysis using ROC method, select randomly 10000 cells for each cell type
Clust_Markers <- FindAllMarkers(rc.integrated,
logfc.threshold=log(2),
min.diff.pct = 0.25,
max.cells.per.ident = 10000,
only.pos=T,
test.use="roc")
# Load in gene ids and names
feature_names <- read_tsv("./supp_data/features.tsv.gz", col_names = c("gene", "Name", "Type")) %>%
select(-Type) %>%
distinct()
# Merge the DE output with corresponding gene id and name
Clust_Markers <- left_join(Clust_Markers, feature_names)
# For each cell type, select the top 3 genes according to AUC score
Clust_Markers %>% group_by(cluster) %>% top_n(3, myAUC)
Clust_Markers$cluster_gene <- paste(Clust_Markers$cluster, Clust_Markers$gene, sep="_")
# store all results
All_clust_markers <- Clust_Markers
# subset based on AUC 0.75
Clust_Markers <- filter(Clust_Markers, myAUC>=0.75)
# Calculate the frequency of genes being detected as marker
(times_per_gene <- Clust_Markers %>%
ungroup() %>%
group_by(gene) %>%
tally())
Clust_Markers <- left_join(Clust_Markers, times_per_gene)
summary(Clust_Markers$n)
# Load libraries for plotting
library(cowplot)
library(ComplexHeatmap)
library(circlize)
library(GeneOverlap)
library(gprofiler2)
library(ggrepel)
library(ggplot2)
options(repr.plot.width = 16, repr.plot.height = 14)
markers_sel <- select(Clust_Markers, gene, cluster)
markers_list <- split(markers_sel, f=markers_sel$cluster)
#this makes list from long df of gene lists - TARGET is what we want to keep
markers_list <- lapply(markers_list, function(x) x[names(x)=="gene"])
# convert each sublist into character and eliminate duplicates
markers_list <- lapply(markers_list, function(x) as.character(unique(x$gene)))
## GeneOverlap
# number of integrated features
genome_size <- 17513L
#compare all lists
gom.self <- newGOM(markers_list, markers_list, genome.size=genome_size)
int <- getNestedList(gom.self, "intersection")
int_matrix <- getMatrix(gom.self, "intersection")
p.val <- getMatrix(gom.self, "pval")
JC <- getMatrix(gom.self, "Jaccard")
# log of p.val for intersection
p.val_log <- -log10(p.val + 1e-200)
olap <- Heatmap(p.val_log,
name = "-log10_pval",
col = colorRamp2(c(0, 200),
c("beige", "red")),
column_title = "Number of shared marker genes",
cluster_rows = T,
cluster_columns = T,
use_raster= FALSE,
show_column_names = TRUE,
show_row_names = TRUE,
show_row_dend = TRUE,
clustering_distance_rows = "pearson",
clustering_distance_columns = "pearson",
show_column_dend = TRUE, cell_fun = function(j, i, x, y, width, height, fill) {grid.text(sprintf("%.0f", int_matrix[i, j]), x, y, gp = gpar(fontsize = 10))
})
# padding - bottom, left, top, right
draw(olap, padding = unit(c(15, 5, 5, 10), "mm"), heatmap_legend_side = "left")
# Load GO terms for Arabidopsis thaliana
cluster_GO <- gost(markers_list, organism = "athaliana", correction_method = "fdr", significant = F, multi_query = F)
cluster_GO_df <- cluster_GO[[1]]
cluster_GO_sig <- filter(cluster_GO_df, p_value<=0.01)
# top terms for each cluster
cluster_GO_sig %>%
filter(source=="GO:BP", intersection_size>=4) %>%
group_by(query) %>%
top_n(5, wt = -p_value) %>%
arrange(desc(p_value)) -> top_GO
GO_n <- cluster_GO_sig %>%
filter(source=="GO:BP", intersection_size>=4) %>%
group_by(term_id) %>%
tally() %>%
arrange(desc(n))
GO_n <- dplyr::rename(GO_n, "n_clusters"=n)
cluster_GO_sig_n <- left_join(cluster_GO_sig, GO_n)
# get all terms for the top ones so that all clusters have values
top_GO_all <- filter(cluster_GO_df, term_id %in% top_GO$term_id)
#spread and plot
top_GO_sel <- select(top_GO_all, query, p_value, term_id, term_name)
spread_GO <- spread(top_GO_sel, key = query, p_value)
spread_GO[is.na(spread_GO)] <- 1
spread_GO_m <- as.matrix(-log10(spread_GO[3:ncol(spread_GO)]))
rownames(spread_GO_m) <- spread_GO$term_name
options(repr.plot.width = 18, repr.plot.height = 20)
GO_hm <- Heatmap(spread_GO_m,
name = "-log10_pval",
heatmap_legend_param = list(title_position="topcenter", color_bar = "continuous"),
col = colorRamp2(c(0, 10),
c("beige", "#e31a1c")),
cluster_rows = T,
cluster_columns = T,
use_raster= FALSE,
show_column_names = TRUE,
show_row_names = TRUE,
show_row_dend = TRUE,
show_column_dend = TRUE,
clustering_distance_rows = "pearson",
clustering_distance_columns = "pearson",
row_names_gp = gpar(fontsize = 12))
# padding - bottom, left, top, right
draw(GO_hm, padding = unit(c(15, 15, 5, 80), "mm"), heatmap_legend_side = "left")
# Load the known markers
(used.markers <- read_csv("./supp_data/Curated_Markers.csv"))
used.markers$cluster_gene <- paste(used.markers$Celltype, used.markers$Locus, sep="_")
# Combine DE output to see whether the known markers are captured
Clust_Markers <- mutate(Clust_Markers, pct.diff=pct.1-pct.2)
Clust_Markers <- arrange(Clust_Markers, desc(pct.diff)) %>%
group_by(cluster) %>%
mutate(pct.diff_rank=dplyr::row_number()) %>%
arrange(desc(avg_diff)) %>%
mutate(avg_diff_rank=dplyr::row_number()) %>%
arrange(desc(myAUC)) %>%
mutate(myAUC_rank=dplyr::row_number()) %>%
mutate(combined_rank_raw=(pct.diff_rank + avg_diff_rank + myAUC_rank)/3) %>%
arrange(combined_rank_raw) %>%
mutate(combined_rank=dplyr::row_number()) %>%
select(-combined_rank_raw) %>%
mutate(known_marker=cluster_gene %in% used.markers$cluster_gene) %>%
arrange(combined_rank)
Clust_Markers
# Store the data
write.csv(Clust_Markers, "./supp_data/Atlas_celltype_markers_ROC.csv", row.names=F)
# Plotting the expression pattern of known markers
used.markers.to.plt <- used.markers
used.markers.to.plt$cluster <- used.markers.to.plt$Celltype
used.markers.to.plt$cluster <- factor(used.markers.to.plt$cluster, levels= c("Quiescent Center", "Stem Cell Niche", "Columella", "Lateral Root Cap", "Atrichoblast", "Trichoblast", "Cortex", "Endodermis", "Pericycle", "Phloem", "Xylem", "Procambium"))
used.markers.to.plt <- arrange(used.markers.to.plt, cluster)
used.markers.to.plt
ct_to_plt_Names <- rev(used.markers.to.plt$Gene)
ct_to_plt_genes <- rev(used.markers.to.plt$Locus)
breaks <- rev(ct_to_plt_genes)
labels <- rev(ct_to_plt_Names)
rc.integrated$celltype.anno <- fct_rev(rc.integrated$celltype.anno)
Idents(rc.integrated) <- "celltype.anno"
options(repr.plot.width=14, repr.plot.height=8)
dot <- DotPlot(object = rc.integrated, features = ct_to_plt_genes, cols = c("white", "red")) + RotatedAxis() + ylab("") +
scale_x_discrete(breaks=breaks, labels=labels)
dot
# top genes to plot
novel_markers <- Clust_Markers %>%
group_by(cluster) %>%
top_n(2, -combined_rank)
novel_markers$cluster <- factor(novel_markers$cluster, levels= c("Quiescent Center", "Stem Cell Niche", "Columella", "Lateral Root Cap", "Atrichoblast", "Trichoblast", "Cortex", "Endodermis", "Pericycle", "Phloem", "Xylem", "Procambium"))
novel_markers <- arrange(novel_markers, cluster)
novel_markers
novel_markers_sel <- select(novel_markers, cluster, gene, Name)
ct_to_plt_Names <- rev(novel_markers_sel$Name)
ct_to_plt_genes <- rev(novel_markers_sel$gene)
(breaks <- rev(ct_to_plt_genes))
(labels <- rev(ct_to_plt_Names))
options(repr.plot.width=10, repr.plot.height=6)
dot <- DotPlot(object = rc.integrated, features = ct_to_plt_genes, cols = c("white", "red")) + RotatedAxis() + ylab("") +
scale_x_discrete(breaks=breaks, labels=labels)
dot
table(rc.integrated$time.celltype.anno)
### combination of time and cell type anno
rc.integrated$cell_time_final <- rc.integrated$time.celltype.anno
rc.integrated$cell_time_final <- gsub("Meristem", "Meristematic", rc.integrated$cell_time_final, ignore.case = FALSE, perl = FALSE,
fixed = T, useBytes = FALSE)
rc.integrated$cell_time_final <- gsub("Elongation", "Elongating", rc.integrated$cell_time_final, ignore.case = FALSE, perl = FALSE,
fixed = T, useBytes = FALSE)
rc.integrated$cell_time_final <- gsub("Maturation", "Mature", rc.integrated$cell_time_final, ignore.case = FALSE, perl = FALSE,
fixed = T, useBytes = FALSE)
# Meristematic-Putative Quiescent Center
rc.integrated$cell_time_final <- gsub("Meristematic-Putative Quiescent Center", "Meristematic-Quiescent Center", rc.integrated$cell_time_final, ignore.case = FALSE, perl = FALSE,
fixed = T, useBytes = FALSE)
table(rc.integrated$cell_time_final)
Idents(rc.integrated) <- "cell_time_final"
DefaultAssay(rc.integrated) <- "integrated"
cell_time_Clust_Markers <- FindAllMarkers(rc.integrated,
logfc.threshold=log(2),
min.diff.pct = 0.25,
max.cells.per.ident = 10000,
only.pos=T,
test.use="roc")
cell_time_Clust_Markers
cell_time_Clust_Markers %>%
group_by(cluster) %>%
tally()
cell_time_Clust_Markers %>%
filter(myAUC>=0.75) %>%
group_by(cluster) %>%
tally()
# store all
all_cell_time_Clust_Markers <- cell_time_Clust_Markers
# keep only markers 0.75 and above
cell_time_Clust_Markers <- filter(cell_time_Clust_Markers, myAUC>=0.75)
(times_per_gene <- cell_time_Clust_Markers %>%
ungroup() %>%
group_by(gene) %>%
tally())
cell_time_Clust_Markers <- left_join(cell_time_Clust_Markers, times_per_gene)
cell_time_Clust_Markers <- left_join(cell_time_Clust_Markers, feature_names)
cell_time_Clust_Markers$cluster_gene <- paste(cell_time_Clust_Markers$cluster, Clust_Markers$gene, sep="_")
cell_time_Clust_Markers <- mutate(cell_time_Clust_Markers, pct.diff=pct.1-pct.2)
cell_time_Clust_Markers
cell_time_Clust_Markers <- arrange(cell_time_Clust_Markers, desc(pct.diff)) %>%
group_by(cluster) %>%
mutate(pct.diff_rank=dplyr::row_number()) %>%
arrange(desc(avg_diff)) %>%
mutate(avg_diff_rank=dplyr::row_number()) %>%
arrange(desc(myAUC)) %>%
mutate(myAUC_rank=dplyr::row_number()) %>%
mutate(combined_rank_raw=(pct.diff_rank + avg_diff_rank + myAUC_rank)/3) %>%
arrange(combined_rank_raw) %>%
mutate(combined_rank=dplyr::row_number()) %>%
select(-combined_rank_raw) %>%
arrange(combined_rank)
cell_time_Clust_Markers
options(repr.plot.width = 16, repr.plot.height = 14)
markers_sel <- filter(cell_time_Clust_Markers, myAUC>=0.75) %>% select(gene, cluster)
markers_list <- split(markers_sel, f=markers_sel$cluster)
#this makes list from long df of gene lists - TARGET is what we want to keep
markers_list <- lapply(markers_list, function(x) x[names(x)=="gene"])
# convert each sublist into character and eliminate duplicates
markers_list <- lapply(markers_list, function(x) as.character(unique(x$gene)))
## GeneOverlap
# number of integrated features
genome_size <- 17513L
#compare all lists
gom.self <- newGOM(markers_list, markers_list, genome.size=genome_size)
int <- getNestedList(gom.self, "intersection")
int_matrix <- getMatrix(gom.self, "intersection")
p.val <- getMatrix(gom.self, "pval")
JC <- getMatrix(gom.self, "Jaccard")
# log of p.val for intersection
p.val_log <- -log10(p.val + 1e-200)
olap <- Heatmap(p.val_log,
name = "-log10_pval",
col = colorRamp2(c(0, 200),
c("beige", "red")),
column_title = "Number of shared marker genes - AUC 0.75+",
cluster_rows = T,
cluster_columns = T,
use_raster= FALSE,
show_column_names = TRUE,
show_row_names = TRUE,
show_row_dend = TRUE,
clustering_distance_rows = "pearson",
clustering_distance_columns = "pearson",
show_column_dend = TRUE, cell_fun = function(j, i, x, y, width, height, fill) {grid.text(sprintf("%.0f", int_matrix[i, j]), x, y, gp = gpar(fontsize = 10))
})
# padding - bottom, left, top, right
draw(olap, padding = unit(c(15, 5, 5, 10), "mm"), heatmap_legend_side = "left")
markers_sel <- filter(cell_time_Clust_Markers, myAUC>=0.75) %>% select(gene, cluster)
markers_list <- split(markers_sel, f=markers_sel$cluster)
#this makes list from long df of gene lists - TARGET is what we want to keep
markers_list <- lapply(markers_list, function(x) x[names(x)=="gene"])
# convert each sublist into character and eliminate duplicates
markers_list <- lapply(markers_list, function(x) as.character(unique(x$gene)))
cluster_GO <- gost(markers_list, organism = "athaliana", correction_method = "fdr", significant = F, multi_query = F)
cluster_GO_df <- cluster_GO[[1]]
cluster_GO_sig <- filter(cluster_GO_df, p_value<=0.01)
# top terms for each cluster
cluster_GO_sig %>%
filter(source=="GO:BP", intersection_size>=4) %>%
group_by(query) %>%
top_n(5, wt = -p_value) %>%
arrange(desc(p_value)) -> top_GO
GO_n <- cluster_GO_sig %>%
filter(source=="GO:BP", intersection_size>=4) %>%
group_by(term_id) %>%
tally() %>%
arrange(desc(n))
GO_n <- dplyr::rename(GO_n, "n_clusters"=n)
cluster_GO_sig_n <- left_join(cluster_GO_sig, GO_n)
# get all terms for the top ones so that all clusters have values
top_GO_all <- filter(cluster_GO_df, term_id %in% top_GO$term_id)
#spread and plot
top_GO_sel <- select(top_GO_all, query, p_value, term_id, term_name)
spread_GO <- spread(top_GO_sel, key = query, p_value)
spread_GO[is.na(spread_GO)] <- 1
spread_GO_m <- as.matrix(-log10(spread_GO[3:ncol(spread_GO)]))
rownames(spread_GO_m) <- spread_GO$term_name
options(repr.plot.width = 18, repr.plot.height = 20)
GO_hm <- Heatmap(spread_GO_m,
name = "-log10_pval",
heatmap_legend_param = list(title_position="topcenter", color_bar = "continuous"),
col = colorRamp2(c(0, 10),
c("beige", "#e31a1c")),
cluster_rows = T,
cluster_columns = T,
use_raster= FALSE,
show_column_names = TRUE,
show_row_names = TRUE,
show_row_dend = TRUE,
show_column_dend = TRUE,
clustering_distance_rows = "pearson",
clustering_distance_columns = "pearson",
row_names_gp = gpar(fontsize = 12))
# padding - bottom, left, top, right
draw(GO_hm, padding = unit(c(15, 15, 5, 80), "mm"), heatmap_legend_side = "left")
cell_time_Clust_Markers$clus_to_split <- cell_time_Clust_Markers$cluster
cell_time_Clust_Markers <- separate(cell_time_Clust_Markers, clus_to_split, into=c("timezone", "celltype"), sep="-")
cell_time_Clust_Markers$cluster <- as.character(cell_time_Clust_Markers$cluster)
cell_time_Clust_Markers$timezone <- factor(cell_time_Clust_Markers$timezone, levels=c("Mature", "Elongating", "Meristematic"))
ct_levels <- names(table(rc.integrated$celltype.anno))
cell_time_Clust_Markers$celltype <- factor(cell_time_Clust_Markers$celltype, levels=ct_levels)
cell_time_Clust_Markers <- arrange(cell_time_Clust_Markers, timezone) %>% arrange(celltype)
cell_time_Clust_Markers
write.csv(cell_time_Clust_Markers, "./supp_data/Atlas_cell_time_combined_ROC_Clust_Markers_20200527.csv", row.names=F)
SCR, MYB36 and CASP1 were selected for the endodermis because they are highlighted later in the manuscript and are highly enrichened. Note that MYB36 is expressed in all developmental stages of endodermis, but is highest in the elongation zone.
For others, I took the top 20 combined rank and gave preference for those that were specific (i.e. n<=1) or where known markers
(Markers_to_plot <- read_csv("./supp_data/Atlas_cell_time_combined_ROC_Clust_Markers_20200527_to_plot.csv") %>% filter(plot=="y"))
(c_levels <- unique(Markers_to_plot$cluster))
ct_to_plt_Names <- Markers_to_plot$Name
ct_to_plt_genes <- Markers_to_plot$gene
table(rc.integrated$cell_time_final)
rc.integrated$cell_time_final <- factor(rc.integrated$cell_time_final, levels=rev(c_levels))
Idents(rc.integrated) <- "cell_time_final"
table(rc.integrated$cell_time_final)
(breaks <- rev(ct_to_plt_genes))
(labels <- rev(ct_to_plt_Names))
options(repr.plot.width=14, repr.plot.height=10)
DotPlot(object = rc.integrated, features = ct_to_plt_genes, cols = c("white", "red")) + RotatedAxis() +
scale_x_discrete(breaks=breaks, labels=labels) + xlab("")
options(repr.plot.width=14, repr.plot.height=10)
dot <- DotPlot(object = rc.integrated, features = ct_to_plt_genes, cols = c("white", "red")) + RotatedAxis() + ylab("") + theme(axis.title.y=element_blank(),
axis.text.y=element_blank(),
axis.ticks.y=element_blank()) +
scale_x_discrete(breaks=breaks, labels=labels) + xlab("")
dot
options(repr.plot.width=14, repr.plot.height=10)
dot <- dot + geom_rect(xmin = 0.5, xmax = 1.5, ymin = 0.5, ymax = 1.5, alpha = 0,fill=alpha("grey",0), color="black") +
geom_rect(xmin = 1.5, xmax = 2.5, ymin = 1.5, ymax = 2.5, alpha = 0,fill=alpha("grey",0), color="black") +
geom_rect(xmin = 2.5, xmax = 5.5, ymin = 2.5, ymax = 5.5, alpha = 0,fill=alpha("grey",0), color="black") +
geom_rect(xmin = 5.5, xmax = 8.5, ymin = 5.5, ymax = 8.5, alpha = 0,fill=alpha("grey",0), color="black") +
geom_rect(xmin = 8.5, xmax = 11.5, ymin = 8.5, ymax = 11.5, alpha = 0,fill=alpha("grey",0), color="black") +
geom_rect(xmin = 11.5, xmax = 14.5, ymin = 11.5, ymax = 14.5, alpha = 0,fill=alpha("grey",0), color="black") +
geom_rect(xmin = 14.5, xmax = 17.5, ymin = 14.5, ymax = 17.5, alpha = 0,fill=alpha("grey",0), color="black") +
geom_rect(xmin = 17.5, xmax = 20.5, ymin = 17.5, ymax = 20.5, alpha = 0,fill=alpha("grey",0), color="black") +
geom_rect(xmin = 20.5, xmax = 23.5, ymin = 20.5, ymax = 23.5, alpha = 0,fill=alpha("grey",0), color="black") +
geom_rect(xmin = 23.5, xmax = 26.5, ymin = 23.5, ymax = 26.5, alpha = 0,fill=alpha("grey",0), color="black") +
geom_rect(xmin = 26.5, xmax = 29.5, ymin = 26.5, ymax = 29.5, alpha = 0,fill=alpha("grey",0), color="black") +
geom_rect(xmin = 29.5, xmax = 32.5, ymin = 29.5, ymax = 32.5, alpha = 0,fill=alpha("grey",0), color="black")
dot
ann_df <- tibble(timezone=as.character(Markers_to_plot$timezone),
celltype= as.character(Markers_to_plot$celltype)) %>%
mutate(row=dplyr::row_number(),
col=rep("A")) %>%
group_by(celltype) %>%
mutate(ct_row=dplyr::row_number()) %>%
ungroup() %>%
mutate(Label=celltype)
ann_df$Label[ann_df$ct_row!=2] <- " "
ann_df$Label[ann_df$celltype=="Stem Cell Niche"] <- "Stem Cell Niche"
ann_df$Label[ann_df$celltype=="Quiescent Center"] <- "Quiescent Center"
ann_df$row <- factor(ann_df$row)
ann_df$timezone <- factor(ann_df$timezone, levels=c("Meristematic", "Elongating", "Mature"))
ann_df
options(repr.plot.width=1, repr.plot.height=10)
(time_ann_p <- ann_df %>% ggplot(aes(x=col, y=rev(row), fill=timezone)) + geom_tile() + scale_fill_manual(values = c("#DCEDC8", "#42B3D5", "#1A237E")) + theme_nothing())
options(repr.plot.width=2, repr.plot.height=10)
ann_df %>% ggplot(aes(x=col, y=rev(row), fill=timezone)) + geom_tile() + scale_fill_manual(values = c("#DCEDC8", "#42B3D5", "#1A237E")) + theme_minimal()
options(repr.plot.width=1, repr.plot.height=10)
ann_df$celltype <- factor(ann_df$celltype, levels = order[sort(match(unique(ann_df$celltype),order))])
color <- palette[sort(match(unique(ann_df$celltype),order))]
(cell_ann_p <- ann_df %>% ggplot(aes(x=col, y=rev(row), fill=celltype)) + geom_tile() + scale_fill_manual(values = color) + theme_nothing())
options(repr.plot.width=1, repr.plot.height=10)
cell_ann_p | time_ann_p
options(repr.plot.width=2.2, repr.plot.height=10)
(text_ann_p <- ann_df %>% ggplot(aes(x=col, y=rev(row), label=Label)) + geom_text(size = 5) + theme(plot.margin = margin(0, 0, 35, 35, "cm")) + theme_nothing())
options(repr.plot.width=4.5, repr.plot.height=10)
(c_ann <- text_ann_p | (cell_ann_p | time_ann_p))
options(repr.plot.width=14, repr.plot.height=10)
combined_marker_plot <- text_ann_p + cell_ann_p + time_ann_p + dot + patchwork::plot_layout(widths = c(0.8,0.2,0.2,4), nrow=1)
combined_marker_plot
options(repr.plot.width=9, repr.plot.height=6)
order <- c("Quiescent Center", "Stem Cell Niche", "Columella", "Lateral Root Cap", "Atrichoblast", "Trichoblast", "Cortex", "Endodermis", "Pericycle", "Phloem", "Xylem", "Procambium", "Unknown")
order <- rev(order)
palette <- c("#BD53FF", "#DED3FE", "#5AB953", "#BFEF45", "#008080", "#21B6A8", "#82B6FF", "#0000FF","#FF9900","#E6194B", "#9A6324", "#FFE119","#EEEEEE")
palette <- rev(palette)
rc.integrated$celltype.anno <- factor(rc.integrated$celltype.anno, levels = order[sort(match(unique(rc.integrated$celltype.anno),order))])
color <- palette[sort(match(unique(rc.integrated$celltype.anno),order))]
cell_dim <- DimPlot(rc.integrated, reduction = "umap", group.by = "celltype.anno", cols=color)
cell_dim[[1]]$layers[[1]]$aes_params$alpha = 1 # can change w/ this if you like. I decided not to so side anno would match
cell_dim
options(repr.plot.width=9, repr.plot.height=6)
time_plt <- DimPlot(rc.integrated,
group.by = "time.anno",
order = c("Meristem","Elongation","Maturation"),
cols = c("#1A237E", "#42B3D5", "#DCEDC8"))
time_plt
options(repr.plot.width=16, repr.plot.height=18)
(cell_dim + time_plt) / combined_marker_plot + patchwork::plot_layout(nrow=2, heights = c(1, 1.4), guides = 'collect')
ggsave("./supp_data/Atlas_cell_dev_marker_dot.pdf", width=16, height=18)
sessionInfo()